Radiation pattern of a classical dipole in a photonic crystal: photon focusing 



Dmitry N. Chigri J* 

Institute of Materials Science and Department of Electrical and Information Engineering, 
University of Wuppertal, Gauss-str. 20, 42091 Wuppertal, Germany 

The asymptotic analysis of the radiation pattern of a classical dipole in a photonic crystal pos- 
sessing an incomplete photonic bandgap is presented. The far- field radiation pattern demonstrates 
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I. INTRODUCTION 

Purcell was the first who pointed out, that the spon- 
taneous emission of an atom or a molecule depends on 
its environment. Since then, an influence of non-trivial 
boundary conditions in the vicinity of an excited atom 
on its emissive properties has been the subject of active 
research 0, 0, ■ Important examples of such an influ- 
ence are an enhancement and inhibition of the sponta- 
neous emission by a resonant environment e.g., mi- 
crocavity. These phenomena were first demonstrated by 
Goy et al. Q and Kleppner 0, respectively, and con- 
tinue to be the subject of intense research not only due 
to their contribution to the better understanding of the 
light matter interaction, but, to a great extend, due to 
the practical importance of controlling the light emission 
process. Light-emitting diodes [EISEl and thresholdless 
lasers [H [H E2 are just a few examples, where the light 
extraction and the spontaneous emission control by mean 
of optical microcavity leads to improved performance. 

Dielectric periodic medium, also called photonic crys- 
tal EE E3, is a good example of non-trivial boundary 
conditions on electromagnetic field. Such an inhomoge- 
neous medium can possess a complete photonic bandgap, 
i.e., a continuous spectral range within which linear prop- 
agation of light is prohibited in all spatial directions. One 
of the consequence is an inhibited spontaneous emission 
for the atomic transition frequency inside the complete 
photonic bandgap 0, 0, UJ. There are no electro- 
magnetic modes available to carry the energy away from 
the atom at complete photonic bandgap frequencies. Al- 
though an existence of complete photonic bandgap usu- 
ally requires a high index (n > 3) dielectric materials 
arranged in a three-dimensional (3D) lattice 0,0], pho- 
tonic crystals are proven to be useful artificial materials 
to modify the light emission even in the absence of com- 
plete photonic bandgap. For example, it was demon- 
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strated that the external quantum efficiency of light- 
emitting diodes can be significantly improved byintro- 
ducing a two-dimensional (2D) photonic crystal |0,E1- 
Another example is a highly directive light source em- 
ploying a 3D photonic crystal [20L l2l| . 

An intrinsic property of photonic crystals is their com- 
plicated photonic band structure, which can be engi- 
neered by choosing an appropriate combination of mate- 
rials and lattice geometry p^.ll4|. Being able to modify 
in purpose the emission rate within a specific spectral 
range and simultaneously in specific directions could add 
a significant flexibility in improving light sources. 

A number of papers were devoted to the study of the 
spontaneous emission in photonic crystals, conside ring 
emission modification using both classica l 1221 1231 124 . l25l 
HI E3 and quantum 0, HI |M US iJIHlsll for- 
malisms. But, to the author knowledge, questions like 
modification of the emission rate in a specific direction 
and modification of the emission pattern due to the pho- 
tonic crystal environment have not been yet addressed. 
Special opportunities in controlling directionality of emis- 
sion exist within spectral ranges of allowed photonic 
bands, where photonic crystals display strong dispersion 
and anisotropy. The consequence of anisotropy is the 
beam steering effect 0,0], which in the essence means, 
that the group velocity direction of the medium's eigen- 
mode does not necessarily coincide with its wave vector 
direction. A beam steering effect known to be a reason 
for the number of anomalies in an electromagnetic beam 
propagation inside a photonic crystal, which are usually 
referred to as superprism or ultrarefractive phenomena 
0, EE HSl- F° r example, an extraordinary large or neg- 
ative beam bending [3g , a beam self-collimation [3?l |3S| 
and the photon focusing [3^. l4fj| were reported. The 
last phenomenon is similar to the phonon focusing, phe- 
nomenon observed in the ballistic transport of phonons 
in crystalline solid [4l| . 

The term phonon focusing refers to the strong 
anisotropy of heat flux in crystalline solid. First observed 
in 1969 by Taylor et al. |42j, phonon focusing is a prop- 
erty of all crystals at low temperatures. The term "fo- 
cusing" does not imply a bending of particle paths, as in 
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the geometrical-optics sense of the term. The physical 
reason of the phonon focusing is the beam steering. In 
particular, waves with quite different wave vectors can 
have nearly the same group velocity, so the energy flux 
associated with those waves bunches along certain crys- 
talline directions. In some special heat flux can 
display intricate focusing caustics, along which flux tends 
to infinity |4l| . This happens when the direction of the 
group velocity is stationary with respect to a small vari- 
ation of the wave vector. 

One can expect, that a similar phenomenon takes place 
in photonic crystals [39l Ecj . An optical cousin of the 
acoustic phenomenon opens a unique opportunity to de- 
sign a caustics pattern on purpose, enhancing and sup- 
pressing emission in specific directions. 

In this paper a description of an angular distribution of 
radiated power of a classical dipole embedded in a pho- 
tonic crystal is presented. It is assumed that only prop- 
agating modes of the photonic crystal contribute to the 
far-field radiation. The emission process is treated using 
an entirely classical model, similar to one in 00. It 
is commonly accepted that a classical description leads 
to the same results as an entirely quantum electrody- 
namical approach |23,|2^|. In the classical description, 
the modification of spontaneous emission is due to the 
radiation reaction of the back-reflected field on the clas- 
sical dipole [43L lil . l45j . Then within the framework of 
the Weisskopf-Wigner approximation 0,^3 , the spon- 
taneous emission rate, T, is related to the classical radi- 
ated power P(r ) = (w/2)Im[d* • E(r )] via T = P/Tluj 
|44j . where d is a real dipole moment, E(ro) is a field in 
the system and ro is the dipole location. A well known 
interpretation of the emission rate modification, as the 
dipole interaction with the out-of-phase part of the radi- 
ation reaction field, follows from that relation [Zsl I45I ] . In 
the classical picture, a non-relativistic Lamb shift is due 
to the dipole interacting with the in-phase part of the 
reaction field liM 1471 and it is also seen to be a purely 
classical effect [4a, • Although, the magnitude of the 
anomalous Lamb shift in a realistic photonic crystal is ac- 
tively and controversially discussed [28|, 123, bMJ, L13 > this 
question is out of the scope of this work. 

A general expressions for the field and emission rate 
of the point dipole radiating in an arbitrary periodic 
medium are reviewed in sections ITU and ITTT1 respectively. 
The evaluation of the asymptotic form of the radiated 
field is given in section llVl In section an angular dis- 
tribution of radiated power is introduced. A modification 
of radiation pattern is discussed in terms of photon fo- 
cusing in section IVII A numerical example of an angular 
distribution of emission power radiated from the point 
isotropic light source is presented in section lYlII for the 
case of a two-dimensional square lattice photonic crys- 
tal of dielectric rods in air. Summary is given in section 



II. NORMAL MODE EXPANSION OF DIPOLE 
FIELD 

In this paper, a general linear, non-magnetic, dielectric 
medium with arbitrary 3D periodic dielectric function, 
e(r) = e(r + R), is studied. Here R is a vector of the 
direct Bravais lattice, R = Y). IjSLj , h is an integer and 
&i is a basis vector of the periodic lattice. It is assumed 
that a medium is infinitely extended in space and that 
no absorption happens. In general, presented approach 
is valid for any inhomogeneous, non-absorbing medium, 
for which a dispersion relations can be found in the form, 
uj = ui(k), numerically or analytically. Here k is the wave 
vector. 

In Gaussian units, Maxwell's equations in such a 
medium have a form 



VXE = -C^' 

VxH=i £ (r)f 

c ot 

V • [e(r)E] - 0, 
V-H = 0. 



47T , 



(1) 

(2) 

(3) 
(4) 



Here, the electric (magnetic) field is denoted by E (H), 
c is a speed of light in vacuum. An electromagnetic field 
is produced by a current source J and the charge den- 
sity is zero, p = 0. Then one can choose the transverse 
(Coulomb) gauge for the vector potential A in the form 



V • [e(r)A] = 0. 



(5) 



The absence of the charge density implies that the scalar 
potential tp is zero. The electric and magnetic fields can 
be written in terms of the vector potential A via: 



E 



I dA 
H = V x A. 



(6) 
(7) 



Combining equations I|6I7I) with Maxwell's equations 
one obtains the wave equation for the vector potential 

A: 

VxVxA + l £ (r)^ = ^J. (8) 

In what follows, a simplest form of the current density J 
is taken: 



J(r, t) = — ioJod5(r — r )e 



(9) 



for a harmonically oscillating dipole with a frequency loq 
and a real dipole moment d, located at the position r 
inside a photonic crystal, switched on at t = 0. 

The field of an arbitrary light source embedded in a 
periodic medium can be constructed by a suitable super- 
position of the medium's eigenwaves: 
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A(v,t)=y2 d 3 k n C nk (t) A„ k (r). (10) 
„ Jbz 

Here A nk (r) and C„ k (i) are the Bloch eigenvector and 
the time-dependent amplitude coefficient of the eigen- 
wave (n, k), respectively. The form of the amplitude co- 
efficient is defined by the particular nature of the light 
source. The integration is performed over the first Bril- 
louin zone (BZ) of the crystal and the summation is car- 
ried out over different photonic bands, where n is the 
band index. 

Eigenwaves A nk (r) satisfy the homogeneous wave 
equation 

V x V x A„ k - 4r £ (r)A„ k = ( n ) 
and also fulfill the orthogonalization, normalization and 



closure conditions given by: 

/ d 3 re(r)A„ k (r)A:, k ,(r) = V6 nn ,5(k - k'), (12) 
Jv 

J d 3 kA Ilk (r)A; ik (r') = l E± S(r - r'), (13) 

where w nk is the Bloch eigenfrequency, V is the volume 
of the unit cell of the crystal, * denotes the complex con- 
jugate and I £± is the identity operator on the subset 
of the e-transverse vector functions as defined in 
The Bloch eigenvector A nk (r) obeys the gauge condition 
V • [e(r)A„ k (r)] — and are therefore transverse with 
respect to this gauge. Equations (|12I13|) ensure that the 
eigenvectors A nk (r) form a complete set of orthonormal 
e-transverse functions. Here any vector that satisfies the 
e-transverse gauge condition 101 is called "e-transverse" 
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The amplitude coefficients C nk (t) can be easily obtained from the wave equation JHJl. Substituting (|10f> into the 
wave equation (JSJ) and using the homogeneous wave equation l|llf) . one obtains 

]T J Bz d 3 k n ( d 2( g (t) + o£ k Cnk(t)) e(r)A nk (r) = W(r,i) 



Then taking the inner product between every term of 
this equation and an eigenvector A n / k /(r), i.e., multiply- 
ing by A*, k ,(r) and integrating over the unit cell of the 
crystal, one finally obtains the differential equation for 
the amplitude coefficients C nk (i) 



d 2 C nk (t) 
dt 2 



- w r 2 lk C„ k (i) 



o 



V 



Ank(ro) ■ d) e 



-lUjQt 



where the orthogonality of the eigenvectors i|12fl and a 
specific form of the source term @ were taken into ac- 
count. Then assuming initial conditions C„ k (0) = 0, one 
has the following solution of this differential equation 



Aircup (A,; k (r ) ■ d) _ 

' v K k - wg) 



(14) 



Finally, the electromagnetic field radiated by a point 
dipole located at r can be represented in terms of Bloch 
normal modes as: 



A(r,i) 



3l , « k (ro)-d) 



d J k 



BZ 



x a nk (r)e lk "-( r ~ ro )e-^ ot , 



(15) 



where the Bloch theorem, A nk (r) = a„ k (r)e lk "' r , have 
been used. 

The integrand in l|15|) has a pole at = , and the 
integral is singular. This is a typical behavior for any 



resonance system, where dissipation is neglected. The 
standard way to regularize the integral is to add a small 
imaginary part to ujq. The result of the integration then 
becomes dependent on the sign of this imaginary part. 
The criterion for determining the sign will be discussed 
below. A regularized integral (|15|l reads 



A(r,t) 



AircujQ \ ^ 



V 

x a„ k (r)e 



BZ 



J3 V KkM ■ d) 



ik„-(r— ro)g— Mat 



(16) 



III. SPONTANEOUS EMISSION RATE 



A light source situated in an inhomogeneous medium 
is immersed in its own electric field emitted at an earlier 
time and reflected from inhomogeneities in the medium. 
By conservation of energy, the decay rate at which en- 
ergy is radiated is equal to the rate at which the charge 
distribution of the source does work on the surrounding 
electromagnetic field. For an arbitrary current density J, 
the radiated power is given by |48| : 



P(t) = - [ d 3 rJ(r,i)-E(r,t), 
Jv 



(17) 



where V is a volume containing a current density source 
J and it is related to spontaneous emission rate via T = 



4 



eZQ nk 



Md 2 k 



,( P V i 

v nk 




(22) 



Now, making use of the time-reversal symmetry of the 
Maxwell's equations, which for a periodic medium im- 
plies that u> n ^ — d> n ,-k |l3j |. one can rewrite the time 
integral in IL'21) to obtain 



dte i(u 2 -ul k )t 



J* - 



7r5(wo-a^ k ). 



FIG. 1: Diagram showing the relations between k-space and 
coordinate space quantities. Iso-frequency contours for fre- 
quencies w„k and u> n k + dui are presented. 



p/huo m 

has 



For the time-averaged radiated power, one 



P = — Re 
2 



/ d 3 rJ*(r,i) -E(r,i) 
Jv 



or, specializing to a point dipole © 



P 



^lm[d* 
2 1 



E(r )] • 



(18) 



(19) 



A well known interpretation of the emission rate modifi- 
cation follows from this relation. Emission rate is mod- 
ified due to the dipole interaction with the out-of-phase 
part of the radiation reaction field 0, 03 ■ 

Consider an excited molecule or atom at a position 
ro in a photonic crystal. Assuming that the presence 
of the molecule does not change the band structure of 
the crystal, the only possible mode it can emit in, is an 
eigenmode of the photonic crystal. In the classical ap- 
proach, the molecule is modeled by a point dipole d JjJJ • 
The radiation reaction field will be given by the normal 
mode expansion 1)16(1 , which is valid for any point r in the 
crystal, which is distinct from (but as close as required 
to) the dipole location r . Such a choice of the radiation 
reaction field corresponds to the Weisskopf-Wigner ap- 
proximation |3lL l46fl . Then the radiated power (emission 
rate) l|19f) of a classical dipole in a photonic crystal is 
given by: 



P = Im 



V 



-£ 



BZ 



3 |A nk (r ) -d 



(20) 



where equation relating E = — (l/c)<9A/9t, was 
used. This integral can be converted to a two- 
dimensional integral over the iso-frequency surface in k- 
space. With the aid of the integral representation 



1 



dre 



X — l JQ 

one can transform the integral i|20|) to the form: 



(21) 



Im 



V ^ 



d 3 k n |A„ k (r ) ■ d| 



BZ 



Then, changing the integration variable to the eigenfre- 
quency cc> nk by use of the relations | V k w„ k | dk — doj„ k 
and d 3 k n = dkd 2 \i n , where d 2 k„ is an element of the iso- 
frequency surface cj„ k = luq (Fig. QJ, equation l|22|) can 
be transformed to: 



V ^ 



A nk (r ) • d 2 2 
awnk — = ; — d(w -w„ k ), 

VkWnk 



where one can carry out the integration over the eigen- 
frequency u nk to obtain finally 



r2, ,2 



V 



£ 



M-dr 



iv„ 



(23) 



where, V nk = V k cj„ k , the group velocity of the eigen- 
wave (n, k) is introduced. 

Formula l)23|) gives the total time-averaged radiated 
power of a dipole situated inside a photonic crystal 
|22t l24j . This result agrees with fully quantum electro- 
dynamical result for spontaneous emission rate of a two- 
level atom within the Weisskopf-Wigner approximation 

m. 



IV. ASYMPTOTIC FORM OF DIPOLE FIELD 

In this section, a radiating dipole field is analyzed in 
the radiation zone. For that, an asymptotic form of the 
integral l|16fl is evaluated and analyzed. In what follows, 
an asymptotic analysis of the Green's function developed 
by Maradudin [4!j for the phonon scattering problem is 
used. 

Using the integral representation l|21|) one can rewrite 
(OH as 



A(r) 



V 



£ 



BZ 



/>oo 

/ dr(a; ik (r ) -d)a„ k (r)e 
Jo 



iF nk (r) 



where 



Ku(t) = K ■ (r - r„) - r(c^ k - w g) 



. (24) 



(25) 



and a limit 7 — > was taken. 

In a typical experiment |x| = |r — r | 3> A, where A is 
the wavelength of the electromagnetic wave. For large |x| 
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an exponential function in the integral l|24[l will oscillate 
very rapidly and one can use the method of stationary 
phase to evaluate the integral. 

The principal contribution to the integral comes from 
the neighborhood of those points in r- and k-space where 
the variation of F n k(r) is the smallest. This means that 
one can set the gradient of the function F n k(r) in k-space 
equal to zero as well as the derivative of the function with 
respect to r. This gives the conditions 

° Fnk 2 2 " (26) 



8t 



0. 



V k F„ k = x - rV k < k = 0. (27) 

Equations (|26I27|) determine the values of r and k„ 
around which the principal contributions to the integral 
(|24|l arise. These points are called stationary points. Fur- 
ther, the stationary points are denoted by r v and k^. 
Assuming that value of the eigenvector a„k (r) is approx- 
imately constant a„k(r) ss a ^k( r ) f° r T close to r„ and 
for the wave vectors close to k^, the integral l(2"4"|l is re- 
duced to the sum of the i nteg rals in the vicinities of the 
stationary points (r^k^) |4g.l50| 

A(r) « ^££(< k (r ).dK k (r) 



V 



d 3 k n / dre tF ^ T \ 



(28) 



Here an extra summation is over all possible solutions of 
Eqs. 1281571 . 

Due to Eq. (126[) the principal contribution to 
the asymptotic behavior of A(r) comes from the iso- 
frequency surface in k-space defined by w 2 k = uiq or 
equivalently defined by w„k = (eigenfrequency w„k 
is positive and real). At the same time, due to Eq. I|27() 
the portion of the iso-frequency surface w„ k = to Q , which 
contributes to the asymptotic field, is the portion near 
the point on this surface where the gradient VkW 2 k is 
parallel to x. One can express the latter condition in an 
alternative fashion. Equation (|27(l can be simplified as: 

x = 2rw„ k V Ilk , 

where V Ilk = Vk^nk is the group velocity of the eigen- 
wave (n, k). So, equation (|27J) just says that the principal 
contribution to the asymptotic behavior of the field A(r) 
at large |x| = |r — ro| 3> A comes from the neighborhood 
of the points k^ on the iso-frequency surface w n k = at 
which the eigenwave group velocity is collinear to obser- 
vation direction x. Since r is positive by definition (|21|l . 
V^ k and x should not only be collinear, but should point 
in the same direction as well, i. e., x • V^ k > 0. 

Assuming that the major contribution comes from the 
regions near the stationary points, one makes a little error 
by extending the integration in (|28|l over all space 

A(r) « l^^^(ar; k (r ).d)a^ k (r) 



V 



d 3 kr, 



dre iF ^\ 



(29) 



Then, the integral over r is simply given by Dirac 8- 
function: 



dre -(^-L) =2 ^ 2 -t4k) 



and one can further convert the volume integration in k- 
space to an integral over the iso-frequency surface w„k = 
u)q. In fact, by using the relations VkW„k| dk = dw n k and 
d 3 k = dkd 2 \t, and integrating over the eigenfrequency 
w n k, the volume integration over k transforms to: 



d 3 k„e ik "( r - r °)% 2 -^ k ) 



<rk n — 

-oo ^0 



ik„-(r— ro) 

|V Jlk | ' 



where V„k = Vk^nk is the group velocity of the eigen- 
wave (n, k). So, the asymptotic form of the field A(r) is 
given finally by: 

A (r) « ^. c ryfe( r »)' d XtW 



V 



|V£ k | 



<b d 2 k„e lk " (r - ro) , 

J ■ — OO 



(30) 



where the comparatively slowly varying function V nk 
was replaced by its value at stationary point k^ and was 
taken outside the integral over k. 

To evaluate the integrals in Eq. 130(1 the analysis of 
the form of the iso-frequency surface in the vicinity of 
one of the stationary points, k^, should be done. It is 
convenient to introduce the local curvilinear coordinates 
£i with the origin at the stationary point and with one of 
the coordinate aligned perpendicular to the iso-frequency 
surface, e.g., £3. One can expand function ^(^1,^2) = 
k„ • x near the stationary point as: 



1 2 



1 2 

« E fivMZk+O (31) 



i,j',fe=l 



where 



/ d 2 h 



ft V =z 



d 3 h 



d&dtjdb 



and x is a unit vector in the observation direction. All 
derivatives are evaluated at the stationary point k^. 

The result of the integration in (|30|l depends on the lo- 
cal topology of the iso-frequency surface near the station- 
ary point. One can generally classify the local topology 
of the surface by its Gaussian curvature. The Gaussian 
curvature K is the product of the two principal curva- 
tures (inverse radii, K\ and Ki) at a point on the sur- 
face, i.e., K — K\Ki. The points on an iso-frequency 
surface can be elliptical, hyperbolic and parabolic. If the 
Gaussian curvature K > 0, the corresponding point on 
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the iso-frequency surface is called elliptical, and if K < 
it is called hyperbolic. For a complex surface, such as 
the iso-frequency surface in figureEJleft, the regions with 
positive and negative Gaussian curvature alternate. The 
surface is parabolic at the borders between regions with 
curvatures of opposite signs, e.g, convex and saddle. The 
lines along which the curvature changes its sign are called 
parabolic lines. The Gaussian curvature at a parabolic 
point is equal to zero. 

Further, the analysis of the asymptotic form of the 
integral (|3(JI) is undertaken, when the stationary points 
are elliptical or hyperbolic. Then in the close vicinity of 
such a stationary point the following expansion holds: 



where only quadratic terms in the expansion (|31(l were 
kept. By choosing the orientation of the local coordi- 
nates £1 and £2 along the main directions of the surface 
curvature at that point k ra = k^ , one can diagonalize the 



(32) 



matrix oHi- . Then 



M£l)&) = ("l^l + a 2&) > «1 = "ll: "2 = "22- 

(33) 

With such a choice of local coordinates in k-space, the 
product = a^a^ determines the Gaussian curvature 
of the iso-frequency surface at the stationary point k„ = 



Using expansion (|33[1 the asymptotic form of the field l|30l) is now given by: 



A(r) 



47TC 



££ 



«*k(ro)-dK k (r) 



ivr ik i 



d£xd& exp ( «Ci 



(34) 



The integral in Eq. (|34|l is calculated simply to be 

d£exp(^£ 2 ) 



/ 2tt 

x \a\ 



exp 



in 



sign (a) 



and an asymptotic form of the vector potential l|15fl at the position r far from the dipole is given by 
A(r, exp (_i (W + J ( S ign K ) + sign(^)))) ^ (A^^o)- d) A- k (r) _ 



8tt 3 

7^ 



(35) 



(36) 



ro 



where A£ k (r) = a^ k (r)e 



ik^ -r 



and summation is over all stat ionary points with x ■ V^ k > 0. 



According to the Eq. I|36(l the electromagnetic field in- 
side photonic crystal represents a superposition of several 
diverging waves, the number of which equals the number 
of stationary phase points on the iso-frequency surface 
^?ik = w o (Fig. 0-left). Each of these waves has its own 
shape and its own propagation velocity. One comment is 
important here, the asymptotic expansion i|36|) describes 
an outgoing wave (k^ ■ x > 0) only if the corresponding 
group velocity is an outward normal to the iso-frequency 
surface u n k — at point k^. It can happen, however, 
that the group velocity becomes an inward normal for 
some frequencies and regions of k-space (Fig. [2J- left). In 
such a case the dot product k^ • x is not positive in the 
asymptotic expansion l|3t)|l and the expansion describes 
incoming waves. In such a situation, one should change 
the sign of the small imaginary part 7 in regularized equa- 
tion O g3: 



A(r) 



Attcloq 

i 

V 



£ 



BZ 



«k(ro) • d) 



x a„ k (r) e * k "-( r - r °) 



(37) 



and proceed as it has been describe above I|24I36[1 . but 
using the integral representation 

1 1 



x + V) 



dre 1 



(38) 



instead of J2U- 



V. ANGULAR DISTRIBUTION OF RADIATED 
POWER 

In this section, the angular dependence of the dipole 
radiated power (|23J) is introduced. 

Using the definition of the solid angle, df2„k = 
d 2 kcosy>/ |k„| 2 , where G?f2„k is the solid angle subtended 
by the surface element d 2 k„, ip is the angle between the 
wave vector k„ and the group velocity V nk = V k t^ n k 
(Fig. nj, on changing the integration variables, one can 
modify equation (|23|l to the form 



P 



£ 



dn 



rik. 



7T 2 tui iA„ k (r ) • dr ik„r 



V 



|v„ k 



cos if 



(39) 
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FIG. 2: Iso-frequency and wave contours. Left: The central 
region of the iso-frequency contour for normalized frequency 
Q = cud/2-KC = d/\ = 0.569 of an infinite square lattice 2D 
photonic crystal made out of dielectric rods placed in vacuum. 
Rods have the optical index 2.9 and radius r = 0.1 5d, where 
d is the period of the lattice (see Sec. IV 111 for details). The 
stationary points, k 1 , k 2 and k 3 , corresponding to the same 
observation direction x are indicated. Right: Corresponding 
wave contour with folds. The shaded and black regions show 
how two equal solid angle sections in coordinate space (right) 
map to widely varying solid angle sections in k-space (left). 
The wave and group velocity vectors with numbers illustrate 
the folds formation of the wave contour. 



where the function enclosed in the brackets defines the 
radiated power of the dipole per solid angle in k-space 



dP 



|A nk (r ) • d| |k ri 



V 



COS if 



(40) 



To derive the angular distribution of radiated power in 
the coordinate space, one should change the integra- 
tion variables in l|39|) from the k-space to the coordinate 
space. 

The k-space distribution of the radiated power (|40|l is a 
function of the k-space direction, given by the polar, #„k, 
and azimuthal, cf> n k, angles of the wave vector k„. The 
direction of energy propagation in a non-absorbing peri- 
odic medium coincides with the group velocity direction 
|5l) . Whereas the coordinate space angular dependence 
of the radiated power is given by the corresponding group 
velocity direction in the coordinate space (9,<fi). Here 6 
and (j) are the polar and azimuthal angles of the group 
velocity in coordinate space. The k-space to the coor- 
dinate space transformation may be expressed formally 
as 



COS0 = /(cOS0„k, 0nk), 
4> = g(cOS0 n k, 4>nk), 



(41) 
(42) 



where the functions / and g are determined from the 
components of the group velocity vector V^ k || x, where 
x is a unit vector in the observation direction. The Ja- 
cobian of the transformation 11411421) 



df dg df dg 



dcost 



d cos 0„k 



(43) 



\ dQ. 
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w nk\ 






-|tfnk|- 1/2 /f 


Id 2 









V 



??k 



dQ nk 



FIG. 3: Diagram to derive the relation between solid angles 
in the k-space and coordinate space. The iso-frequency con- 
tour for frequency uj nk is presented. The Jacobian of the 
transformation I|41I42|I is given by the ratio dtt/dQ n k. By 
the definition of the solid angle, the solid angle in k-space is 
d£l n k = d kcos</?/ |k„j 2 , while the corresponding solid angle 
in coordinate space is dtl = d 2 k\K nk \. That gives the Jaco- 
bian Jnk = |k„| 2 l-fifnkl / cos tp. Here tp is an angle between the 
wave vector and the group velocity vector. d 2 k is the surface 
element of the iso-frequency surface. 



relates a small solid angle in the coordinate space with 
the corresponding solid angle in k-space via 



dQ = d(cos 



Jnkd(cos0„ k )#nk = Jnh.dfl n k. (44) 



According to the results presented in the Section llVl dif- 
ferent wave vectors can result in the group velocity with 
same direction in coordinate space. That means that the 
following equation 



1 

J nk 



-dQ 



should hold for each stationary wave vector, which sat- 
isfies x • VJ^ k > 0. Changing the integration variables in 
(|39jl one should then sum individual contributions from 
all these wave vectors: 



P = EE 

n v 



4 71 



dQ 



A 2 |A^ k (r )-dn k-| 
V 



Jnk l^nkl 



COS if 



(45) 

The geometrical relationship between solid angles in 
k-space and coordinate space (Fig. |3J) results in the fol- 
lowing formula for the Jacobian (|43[1 



u til 



\Kt, 



, lk | / cos Lp . 

Then, equation 145(1 can be transformed to the form: 



P 



* EE 



^ 2 |A- k (r )-d| 
V IV^kll^kl 



(46) 



where V^ k = VkW„ k is the group velocity and deter- 
mines the Gaussian curvature of the iso-frequency surface 
at the stationary point k„ = k^. Finally, the radiated 
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power of the dipole per solid angle in coordinate space is 
given by the function enclosed in the brackets in i|46|) 



dP 

In 



EE' 



K k (r )-d| 



V V 



nk 



K 



(47) 



Formula (|47|) provides a simple route to calculate an 
angular distribution of radiated power of the point dipole 
© inside a photonic crystal. It can be interpreted as a 
decay rate, at which the dipole transfers energy to the 
electromagnetic waves with the group velocity in the ob- 
servation direction. Then, (dT/dft) — (dP/d£l)/KuiQ is 
related to the probability of the radiative transition of 
an excited atom with emitting a photon traveling in the 
given observation direction. 

Basically, formulae I|46I47|I involve calculations of the 
Bloch wave vectors k^, ending at the iso- frequency sur- 
face w„k = wo, the corresponding group velocity vectors 
V^ k , the Gaussian curvature of the iso- frequency surface 
^nk an d the local coupling strength of the dipole mo- 
ment with a Bloch eigenwave (n,k), given by the factor 
A^ k (ro) • d|. The primary difficulty in obtaining the co- 
ordinate space distribution of radiated power (dP/dtt) 
(|47|l is that the wave vector, the group velocity and the 
Gaussian curvature are all functions of the k-space di- 
rection. Whereas an angular dependence of the radiative 
power (dP/dfl) is given by the corresponding group ve- 
locity direction (6,<fi). To calculate the radiated power 
(dP/dil) (14TB one should take an inverse of the map- 
ping H41I42|I . This inverse is not necessarily unique. In 
the case of multiple stationary points il26l27[) . one direc- 
tion (6, </>) results from several different k-space directions 
($k,0k) (Fig- 13- This requires that the inversion of the 
mapping l|4H42[l must be done point-by-point. 

As a simple exercise, formula l|47f) is applied here to 
calculate an angular distribution of power radiated by a 
dipole in free space. The wave vector and the group ve- 
locity in free space are parallel and their values are simply 
given by |k| = ujq/c and c, respectively. The Gaussian 
curvature of the iso-frequency surface is a square of the 
inverse wave vector 1/ |k| 2 . And the appropriate normal 
modes are plane waves 



A„ k (r) 



V 



(2tt) 3 



<l?ik , 



where a„k is a polarization vector orthogonal to the wave 
vector k. Then, the radiated power is given by (|47|l 



8^' d '- 2 < 



(48) 



'dP\ 
dfl J f 

/ free 

yielding the usual results for radiation pattern in free 
space |48j . 

VI. PHOTON FOCUSING 

The factor |A^ k (r ) • d| 2 in relation l|47|l . giving the 
coupling strength of dipole moment with the photonic 



crystal eigenmode at the dipole position, can display a 
complex angular behavior, which depends on eigenmode 
structure and dipole orientation with respect to the crys- 
tal lattice. To study the net result of the influence of the 
photonic crystal on the radiation pattern of the emitter, 
it is instructive to model an isotropic light source pro- 
ducing a uniform distribution of wave vectors. Moreover, 
an isotropic point source is usually a good model for a 
common experimental situation of emitters with random 
distribution of dipole moment (dye molecules |53 . l53ll54| . 
quantum dots H^Eil, etc.). Then, the radiated power 
(|47|l should be averaged over the dipole moment orienta- 
tion, which simply yields a factor of |d| 2 /3 



dP 

~dh~ 



EE 



(27TC) 3 |A- k (r )| 2 



(49) 



Here the result was normalized to the radiated power 
in free space. Now, the factor |A^ k (r )| 2 gives a field 
strength at the source position and has no angular de- 
pendence. So, the radiation pattern of a point isotropic 
emitter is defined by 



§,) -EE re 



(50) 



The radiated power (|50|l is proportional to the inverse 
group velocity, |V^ k | , and to the inverse Gaussian cur- 
vature, |-?^ k | -1 of iso-frequency surface. A large en- 
hancement of emission rate is expected when the group 
velocity is small. This can be interpreted as a conse- 
quence of the long interaction time of the emitter and 
the radiation field j5|| IBH l58| . In a similar fashion, a 
small Gaussian curvature formally implies an enhance- 
ment of radiated power along a certain observation di- 
rection. While spontaneous emission enhancement due 
to a small group velocity involves non-linear interaction 
of radiation and emitter, the enhancement due to a small 
Gaussian curvature is a linear phenomenon related to the 
anisotropy of the photonic crystal and is a result of the 
beam steering effect. Being a measure of the rate, with 
which emitter transfers energy in photons with a given 
group velocity, radiated power l|50(l will be enhanced if 
many photons with different wave vectors reach the same 
detector. The enhancement of the radiated power, which 
is due to the small Gaussian curvature is called 'photon 
focusing |39l |40| and has a major influence on radiation 
pattern of a point source in a photonic crystal. 

The physical picture of the photon focusing can 
be illustrated in the following manner (Fig. [2J. An 
iso-frequency surface of an isotropic and homogeneous 
medium is a sphere. There is only one stationary point 
with x • V^ k > and thus only one wave propagating 
in the given direction. Figure 0-left is an example of a 
part of the actual iso-frequency contour of a 2D photonic 
crystal made out of dielectric rods placed in vacuum (see 
Sec. IVIll for further details). The anisotropy of the crys- 
tal implies a complex non-spherical iso-frequency surface, 
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which can have several stationary points with x • V^ k > 
(Fig. left). Several waves can propagate in a given di- 
rection inside a photonic crystal. It is illustrative to con- 
struct the wave surface in coordinate space. To construct 
the wave surface one should plot a ray in the observation 
direction x starting from the point source position and 
having the length of the group velocity |V^ k |. An ex- 
ample of the wave contour is presented in figure fright. 
The existence of multiple stationary point implies that 
the wave surface is a complex multivalued surface pa- 
rameterized by wave vector k„. Figure [21 illustrates how 
this can result in a fold of the wave surface. 

In the vicinity of the parabolic point with zero Gaus- 
sian curvature an iso- frequency surface is flat. That im- 
plies, that a very large number of eigenwaves with wave 
vectors in the vicinity of a parabolic point have nearly 
the same group velocity, contributing to the energy flux 
in the direction parallel to that group velocity. In the 
figure |2 it is illustrated by mapping two equal solid an- 
gle sections along different observation direction in the 
coordinate space onto the corresponding solid angle sec- 
tions in k-space [5^. The black solid angle section in 
coordinate space maps onto a single smaller solid angle 
section in k-space implying a "defocusing" of the energy 
flux. The shaded solid angle section in coordinate space, 
which crosses three different branches of the wave con- 
tour, maps onto two different and larger solid angle sec- 
tions in k-space implying enhancement ("focusing") of 
the energy flux in this group velocity direction. This re- 
sults in strongly varying angular distribution of emission 
intensity with sharp singularities (caustics). 



VII. NUMERICAL EXAMPLE: 
CRYSTAL 



2D PHOTONIC 



In this section the theoretical approach developed in 
the previous sections is applied to the numerical calcu- 
lation of the radiation pattern of a point source placed 
inside a 2D photonic crystal. A point source is situated 
inside the crystal and it produces an isotropic and uni- 
form distribution of wave vectors k n with the frequency 

UJ . 

An infinite 2D square lattice of dielectric rods in vac- 
uum (Fig. 0J is considered in the case of in-plane propa- 
gation. Consequently, the problem of an electromagnetic 
wave interaction with a 2D photonic crystal is reduced 
to two independent problems, which are called TE and 
TM, when the electric or magnetic field is parallel to the 
axis of the rods. In the illustrative example presented 
in this section, all numerical calculations have been per- 
formed for TM modes of the crystal. The photonic band 
structure of the crystal made of the rods with the re- 
fractive index n = 2.9 is presented in the figure The 
band structure has been calculated using the plane wave 
expansion method [6(i|. 

In the figure|S]iso-frequency contours of the crystal are 
presented for two frequencies belonging to the first pho- 




FIG. 4: Photonic band structure of TM modes for the square 
lattice photonic crystal with refractive index of the rods n = 
2.9, the lattice constant d and radius of the rods 0.15d. The 
frequency is normalized to f2 = uid/2nc = d/X. c is the speed 
of light in vacuum. Insets show the first Brillouin zone of the 
crystal with the irreducible zone shaded light gray (left) and 
a part of the lattice (right). 




FIG. 5: Iso-frequency contours of the square lattice photonic 
crystal for the normalized frequencies fl = 0.31 (dashed line) 
and SI = 0.34 (solid line). Parabolic point are marked by the 
black dots. The first Brillouin zone of the lattice is plotted 
in order to show the spatial relation between zone boundary 
and iso-frequency contours. 



tonic band (Fig. 0}. To plot an iso-frequency contour, 
the photonic band structure for all wave vectors within 
the irreducible BZ was calculated and then the equation 
w(k) = loq was solved for a given frequency too. Fre- 
quencies have been chosen below (fi = 0.31) and above 
(ft = 0.34) the low edge frequency of the stopband in 
the rX direction of the crystal. The iso-frequency con- 
tours below and above the stopband edge frequency show 
an important difference. As the frequency stays below 
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the stopband, an iso-frequency contour is closed and al- 
most circular (Fig. |SJ . The corresponding wave contour 
(see Section IvTl for definition) is presented in the figure 
|SJ To calculate the group velocity, the plane wave ex- 
pansion method |60| and the Hellmann-Feyman theorem 
were used. The group velocity |V£ k | and the Gaussian 
curvature |-K"^ k | of the iso-frequency contours are rela- 
tively slow function of the wave vector. The Gaussian 
curvature does not vanish for any wave vector. This im- 
plies a small anisotropy in the energy flux inside the crys- 
tal. 

To find how a radiated power varies in coordinate 
space, one should calculate the group velocity and the 
Gaussian curvature on the iso-frequency contour w(k) = 
ujq as functions of an angle in coordinate space. As the 
wave contour is single valued function, the inverse of the 
mapping H41I42J1 from k-space to coordinate space is one- 
to-one and can be easily done. In the figure [3 the polar 
plot of radiated power is presented, which shows small 
amount of anisotropy. The angular distribution of the 
radiated power possesses four-fold rotational symmetry 
of the crystal. 

With increase of the frequency up to the stopband, the 
topology of the iso-frequency contour abruptly changes. 
The stopband developed in the TX direction and the 
iso-frequency contour becomes open (Fig. [5j|. This 
topology changes result in complex contour with al- 
ternating regions of different Gaussian curvature sign. 
Parabolic points, where the Gaussian curvature vanishes, 
are marked by black dots in the figure As it has been 
discussed in section I VII vanishing curvature results in 
the folds of the wave contour. The wave contour cor- 
responding to the iso-frequency Q — 0.34 is presented 
in the figure |S| A pair of the parabolic points in the 
first quarter of the Brillouin zone results in a cuspidal 
structure of the wave contours in the first quarter of the 
coordinate space. This dramatically increases anisotropy 
of the energy flux. 

The folds in the wave contours yields that inverse of 
the mapping I|41I42[1 from k-space to coordinate space 
is not one-to-one any more. To apply the formula H50|) 
to calculate an angular distribution of radiated power 
in such a case, one should proceed as follows. At first, 
the Gaussian curvature as a function of the wave vector 
should be calculated. Then, wave vectors and group ve- 
locities corresponding to the parabolic points on the iso- 
frequency surface should be found. An inversion of the 
mapping (|41I42() should be calculated separately for each 
of the branches of the wave contour. The total radiated 
power is a sum of the different contributions from these 
branches. In the figure|§]the polar plot of radiated power 
(|50|l corresponding to the normalized frequency £1 = 0.34 
is presented. The energy flux is strongly anisotropic for 
this frequency, showing relatively small intensity in the 
directions of the stopband, and infinite intensity (caus- 
tics) in the directions of the folds. 

To substantiate this behavior, finite difference time do- 
main (FDTD) calculations were done [6lll62T|. The sim- 




-i.o p I 

-i.o • •-• • (M) 11 M) 



FIG. 6: Wave contour corresponding to the normalized fre- 
quency Q = 0.31. The group velocity is plotted in the units 
of the speed of light in vacuum. High symmetry directions of 
the square lattice are specified. 
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FIG. 7: Angular distribution of radiative power correspond- 
ing to the normalized frequency £1 — 0.31. High symmetry 
directions of the square lattice are specified. 



ulated structure was a 50 x 50 lattice of dielectric rods 
in vacuum. The crystal is surrounded by an extra 5d 
wide layer of a free space. The simulation domain was 
discretized into squares with a side A = d/32. The total 
simulation region was 1920 x 1920 cells plus 8-cell wide 
perfectly matched layer (PML) (6j|. The poin t isotropic 
light source was modeled by a soft source |6l], |63| with a 
homogeneous spacial dependence and sinusoidal tempo- 
ral dependence of the signal. All FDTD calculations was 
performed with a commercial tool |64| . 

In figure 1101 the map of the modulus of the Poynting 
vector field is shown. The point source is placed in the 
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FIG. 8: Wave contour corresponding to the normalized fre- 
quency Q = 0.34. The group velocity is plotted in the units 
of the speed of light in vacuum. The directions corresponding 
to the folds of the wave contour are shown. 




FIG. 9: Angular distribution of radiative power corresponding 
to the normalized frequency Q = 0.34. The directions of 
infinite radiative power (caustic) coincide with the directions 
of the folds of the wave contour (Fig.|HJl. 



middle of the crystal. The field map is shown for one 
instant time step. The snap-shots were captured after 
10000 time steps, where the time step was 4.38 x 10~ 17 s 
(0.99 of the Courant value). The structure of the crystal 
is superimposed on the field map. From figure lTUI one can 
see, that the emitted light is focused in the directions co- 
inciding with the predicted directions of the folds (black 
lines). 

In figures 1111131 a more complicated example of the 
anisotropy of a photonic crystal is presented. Iso- 
frequency contours for three frequencies belonging to the 




I ■ 

0.0 0.0125 0.25 

FIG. 10: FDTD calculation. Map of the modulus of the 
Poynting vector field for a 50 x 50 rod photonic crystal excited 
by a point isotropic source with the normalized frequency 
Q = 0.34. The location of the crystal in the simulation is 
shown together with asymptotic directions of photon focus- 
ing caustics. 



second photonic band of the crystal are plotted in the fig- 
ure ^2 While iso- frequency contours for the normalized 
frequencies = 0.55 and Q = 0.58 have non-vanishing 
Gaussian curvature for all wave vectors leading to only 
limited anisotropy of the energy flux, the iso-frequcncy 
contour for the normalized frequencies Q = 0.565 displays 
several parabolic points. Moreover, the iso- frequency 
contour consists of two branched with slightly different 
shapes (solid and dashed lines in the figure ITT|l . Two 
branches yield two wave contours with cuspidal folds in 
coordinate space (Fig. I12|l . Applying the formula (|50|l to 
the radiative power calculation, one should sum over con- 
tributions coming from all branches of the wave contours 
in coordinate space. An angular distribution of radiative 
power for the normalized frequencies f2 — 0.565 is pre- 
sented in the figure^H Within the first quarter of the co- 
ordinate space, four caustics with infinite radiated power 
present in the energy flux corresponded to four parabolic 
points on two branches of the iso-frequency contours. 



VIII. SUMMARY 

In this paper, by analyzing a dipole field in the ra- 
diation zone it was shown, that the principal contribu- 
tion to the far-field of the dipole radiating in a photonic 
crystal comes from the regions of the iso-frequency sur- 
face in the wave vector space, at which the eigenwave 
group velocity is parallel to observation direction x. It 
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FIG. 11: Iso-frequency contours of the square lattice photonic 
crystal for the normalized frequencies Q = 0.55 (dotted line), 
fl — 0.565 (solid and dashed line) and £1 = 0.58 (dashed- 
dotted line). Two branches of the iso-frequency contour of 
Q — 0.565 is plotted as solid and dashed lines. Parabolic 
point are marked by the black dots. The first Brillouin zone 
of the lattice is plotted in order to show the spatial relation 
between zone boundary and iso-frequency contours. 



was also shown that anisotropy of a photonic crystal re- 



veals itself in the strongly non-spherical wave front lead- 
ing to modifications of both far-field radiation pattern 
and spontaneous emission rate. By systematic analysis 
of the Maxwell equations a simple formula to calculate 
an angular distribution of radiated power due to a point 
dipole placed in a photonic crystal was derived. The for- 
mula only involves calculations of the wave vectors, the 
group velocity, the coupling strength of the dipole mo- 
ment with the field and the Gaussian curvature on the 
iso-frequency surface corresponding to the frequency of 
the oscillating dipole. That can be done by simple plane 
wave expansion method and is not computationally de- 
manding. A numerical example was given for a square- 
lattice 2D photonic crystal. It was shown by applying 
developed formalism and substantiated by FDTD calcu- 
lations, that if a dipole frequency is within a partial pho- 
tonic bandgap, a far-field radiation pattern is strongly 
modified with respect to the dipole radiation pattern in 
vacuum, demonstrating suppression in the directions of 
the spatial stopband and enhancement in the direction 
of the group velocity, which is stationary with respect to 
a small variation of the wave vector. 
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